Initial Questions to Explore

Is there a correlation of number of hunters to number of elk in each unit? If not, I would certainly be interested in which units have the highest ratio of Elk to Hunters

NOTICE that I am only looking at the general rifle hunting seasons on public land. There are also hunters in Archery, Muzzleloader, Private Land, Ranching for Wildlife, etc.

Setup

setwd("~/_code/colorado-dow/Phase I - Descriptive Analytics")

Load required libraries for wrangling data, charting, and mapping

library(plyr,quietly = T)
library(dplyr,quietly = T)
library(ggplot2, quietly = T)

Set our preferred charting theme

theme_set(theme_minimal())

Run script to get hunter data

source('~/_code/colorado-dow/datasets/read colorado dow pdf.R', echo=F)

Table of the data

COElkRifleAll

Run script to get elk population data

source('~/_code/colorado-dow/datasets/read colorado dow population estimates.R', echo=F)

Table of the data

COElkPopulationAll

Join, fill blanks with NA

COElkHuntersPopulation <- full_join(COElkRifleAll,COElkPopulationAll)
## Joining, by = c("Unit", "Year")

Remove years with no data

COElkHuntersPopulation <- filter(COElkHuntersPopulation, !is.na(Population) & !is.na(Hunters))

Statewide Elk Hunters and Elk Population

First lets look at the entire state as a whole

COElk.Hunters.Population.Statewide <- summarise(group_by(COElkHuntersPopulation,Year,Unit),
                                   Hunters = sum(Hunters,na.rm = T),
                                   Population = sum(Unit_Pop,na.rm = T))

COElk.Hunters.Population.Statewide <- summarise(group_by(COElkHuntersPopulation,Year),
                                   Hunters = sum(Hunters,na.rm = T),
                                   Population = sum(Population,na.rm = T))

# ignore years with missing data
COElk.Hunters.Population.Statewide <- filter(COElk.Hunters.Population.Statewide, Hunters > 0 & Population > 0)

ggplot(COElk.Hunters.Population.Statewide, aes(Hunters,Population,label=Year)) +
  geom_point(size=1) +
  geom_text(hjust=-.2,size=4) +
  geom_smooth(method = 'lm',se=F) +
  labs(title="Statewide Elk Hunters vs Elk Population",subtitle="by Year", caption="source: cpw.state.co.us")

Statewide there is a general relationship to the number of hunters per Elk population. This would make sense as one of the roles of CPW is to regulate the number of hunters to maintain some target Elk Population.

Elk per Hunter by Unit

I’d like to know the distribution of Elk per Hunter across the state.

run script to get unit boundaries so we can draw them on a map

source('~/_code/colorado-dow/datasets/coordinate locations of cpw hunt units.R', echo=F)
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/psarnow/_code/colorado-dow/datasets/CPW_GMUBoundaries/BigGameGMUBoundaries03172015.shp", layer: "BigGameGMUBoundaries03172015"
## with 185 features
## It has 12 fields
## Integer64 fields read as strings:  GMUID

Get a statemap with some roads on it

roaddata <- rgdal::readOGR("~/_code/colorado-dow/datasets/ne_10m_roads/ne_10m_roads.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/psarnow/_code/colorado-dow/datasets/ne_10m_roads/ne_10m_roads.shp", layer: "ne_10m_roads"
## with 56601 features
## It has 29 fields
## Integer64 fields read as strings:  scalerank question
USAroads <- roaddata %>% subset(.,sov_a3 == "USA" & type == "Major Highway")
# I need to convert to data frames so that I can use the data with ggplot2.
USAroads <- fortify(USAroads)
Unitboundaries <- shapefile %>% fortify(region = "GMUID")

Unitboundaries2 <- merge(Unitboundaries, shapefile@data, by.x = 'id', by.y = 'GMUID')
Unitboundaries2$Unit <- as.character(Unitboundaries2$id)

# get min/max of long/lat for zooming
longset <- c(min(Unitboundaries2$long),max(Unitboundaries2$long))
latset <- c(min(Unitboundaries2$lat),max(Unitboundaries2$lat))
COroads <- filter(USAroads, long > longset[1] & long < longset[2])
COroads <- filter(COroads, lat > latset[1] & lat < latset[2])

COElk.Hunters.Population.Unit <- summarise(group_by(COElkHuntersPopulation,Year,Unit),
                                           Hunters = sum(Hunters,na.rm = T),
                                           Population = sum(Unit_Pop,na.rm = T),
                                           Elk_per_Hunter = Population / Hunters)

Year2016 <- filter(COElk.Hunters.Population.Unit, Year == "2016")
ElksperHunterstoPlot <- left_join(Unitboundaries2,Year2016, by=c("Unit"))
ggplot(ElksperHunterstoPlot, aes(long, lat, group = group)) + 
  geom_polygon(aes(fill = Elk_per_Hunter),colour = "grey50", size = .2) + #Unit boundaries
  geom_path(data = COroads,aes(x = long, y = lat, group = group), color="#3878C7",size=2) + #Roads
  geom_text(data=data_centroids,aes(x=longitude,y=latitude,label = Unit),size=3) + #Unit labels
  scale_fill_distiller(palette = "Purples",direction = 1,na.value = 'grey') +
  xlab("") + 
  ylab("") +
  theme(panel.background = element_rect(fill='white')) +
  theme(panel.grid.major= element_blank()) +
  theme(panel.grid.minor= element_blank()) +
  labs(title="2016 Colorado Elk per Hunter", caption="source: cpw.state.co.us")

Conclusion

There are about a dozen units that have much higher ratios of Elk to Hunter. They could be edge cases though. Maybe they are apart of a herd DAU but most of the herd resides in a neighboring hunt unit. Maybe these ratios are legit, and they are hard to get a license in (requires many preference points). Maybe many people apply to hunt there, but the draw is limited and there is a low chance of success (requires preference points). Should investigate these top ranked units.

FUTURE Phase II Diagnostic question – Why do Units 140,20,10,8,471,691,391,851,40 have higher Elk to Hunter ratios than other units?

After looking at this data, I’m also thinking that I’m probably more interested in the success of the hunters regardless of how many elk or how many hunters are in the units. Let’s investigate hunt success next.